Spatiotemporal monitoring of the rare northern dragonhead (Dracocephalum ruyschiana, Lamiaceae) — SNP genotyping and environmental niche modeling herbarium specimens

Abstract The species we have studied the spatiotemporal genetic change in the northern dragonhead, a plant species that has experienced a drastic population decline and habitat loss in Europe. We have added a temporal perspective to the monitoring of northern dragonhead in Norway by genotyping herbarium specimens up to 200 years old. We have also assessed whether northern dragonhead has achieved its potential distribution in Norway. To obtain the genotype data from 130 herbarium specimens collected from 1820 to 2008, mainly from Norway (83) but also beyond (47), we applied a microfluidic array consisting of 96 SNP markers. To assess temporal genetic change, we compared our new genotype data with existing data from modern samples. We used sample metadata and observational records to model the species' environmental niche and potential distribution in Norway. Our results show that the SNP array successfully genotyped all included herbarium specimens. Hence, with the appropriate design procedures, the SNP array technology appears highly promising for genotyping old herbarium specimens. The captured genetic diversity correlates negatively with distance from Norway. The historical‐modern comparisons reveal similar genetic structure and diversity across space and limited genetic change through time in Norway, providing no signs of any regional bottleneck (i.e., spatiotemporal stasis). The regional areas in Norway have remained genetically divergent, however, both from each other and more so from populations outside of Norway, rendering continued protection of the species in Norway relevant. The ENM results suggest that northern dragonhead has not fully achieved its potential distribution in Norway and corroborate that the species is anchored in warmer and drier habitats.


| INTRODUC TI ON
Loss of biodiversity is one of the great challenges facing our society today. In order to predict how climatic and other changes will affect biodiversity, holistic knowledge is needed, from the level of genes to ecosystems. Drivers of biodiversity act both across space and through time. Holistic biodiversity studies should therefore include both spatial and temporal data of various kinds. Contemporary genetic, distributional, and ecological data provide a spatial snapshot across the latest generations only. In the absence of real-time historical genetic data, various methods have been developed for approximating the past and predicting the future impacts of change (i.e., the temporal aspect) based on contemporary data alone. Liberating genetic, distributional, and ecological data from archived biological specimens, however, would enable to create a solid base to study temporal biodiversity dynamics in real time.
Archived biological collections, such as herbaria, fungaria, and seed, culture, in vitro, tissue, and DNA collections, contain expertcurated specimens and associated information (i.e., metadata) collected throughout the world, some of which are several 100 years old. Such scientific collections provide verifiable records of the existence of an organism at a given time and place. With the exception of a few studies, such biological specimen archives have remained a largely untapped resource to study trajectories and trends of genetic diversity (Andrew et al., 2018;Bieker & Martin, 2018), one of the main reasons being recalcitrant DNA in old specimens. Post-mortem degradation of DNA is an inherent trait and unending process of biological materials challenging the usability of archived biological specimens in DNA studies (Allentoft et al., 2012).
With new developments in genomic approaches, genetic data can now be liberated from historical specimens enabling their use as a source for genome-scale biodiversity studies (museomics; Besnard et al., 2018). A second challenge is the lack of standardized, cost-and time-efficient methods for capturing genomic data from historical DNA. Most of the genome-scale approaches currently in use (e.g., shotgun deep sequencing, genome skimming, targeted DNA capture, and de novo organellar genome assembly; see Kistler et al., 2020, and references therein) are still both laborious and expensive and thus are, at this time, of less practical use in large-scale species monitoring and biodiversity assessments. As biodiversity dynamics have no borders or fixed scales, informative biodiversity research calls for large-scale assessments. Hence, there is a need for cheaper and more effective solutions for large-scale biodiversity assessments.
Combining museomics (i.e., genomic analysis of museum specimens) with eco-informatics (i.e., analyses of specimen occurrence data and ecological information) promises to be a fruitful integration of scientific domains, enabling more holistic species knowledge to guide monitoring efforts, in which knowledge about a species' behavior in relation to external forces over time is key. While genomic data can provide an evolutionary framework and delimit units at which selection is operating, georeferenced herbarium records provide basic occurrence data that can be used to understand, predict, and map species distributions and examine past phenological trends and even species interactions (e.g., Meineke et al., 2018).
Occurrence data combined with biotic and abiotic data may further reveal key predictors of the species distributions by surveying various potentially explanatory variables. At times when specimen label information was available only from local, internal databases and often required assistance from collection staff, few researchers made the effort to collect such data. With recent global digitization initiatives, distributional and ecological specimen label information are rapidly becoming readily available through public repositories (e.g., GBIF.org, 2020). Such evolution-ecology integration for species monitoring has been practised for some time (e.g., Bendiksby et al., 2014;Carlsen et al., 2012;Nygaard et al., 2021). Herbarium specimens were an essential data source in these studies, but the temporal dimension was not specifically addressed.
Plants are key components of biodiversity, contributing to ecosystem resilience and services that we depend upon. The world's 178 herbaria (archived collections of plants) contain about 390 million specimens collected throughout the world for more than 350 years (Thiers, 2020). In the present study, we add a historic level to the species monitoring of the flowering plant species, Dracocephalum ruyschiana L. (northern dragonhead; Lamiaceae) by testing a microfluidic-based SNP genotyping array on old herbarium specimens. This approach has recently been applied by others to historical materials of, for example, tropical tree species and salmon (Finch et al., 2020;Johnston et al., 2013;Östergren et al., 2021).
Northern dragonhead is considered a remnant of the glacial steppe flora in Europe with its westernmost occurrences in France and Scandinavia (Lazarević et al., 2009). The typical northern dragonhead habitats throughout its distribution have been, and still are, subject to alteration and destruction by for instance succession (due to ceases of traditional agricultural use like grazing and mowing), increased feralization (due to intensification of farming), and development of infrastructure (due to e.g., use of herbicides along train rails; Norwegian Directorate for Nature Management, 2010;Økokrim, 2013). As future conflicts seem inevitable, it was decided that all conservation options should be assessed, including ex situ conservation (IUCN/SSC, 2014). Translocation is one such method, which has been successfully performed on an entire northern drag- More than 25% of the total European population of northern dragonhead is found in Norway, where it is one of three vascular plants that have a separate percept of law with an action plan for conservation (Lovdata, 2011;Norwegian Directorate for Nature Management, 2010). Northern dragonhead is listed as vulnerable on the Norwegian Red List of 2021 with an estimated 20-40% of populations having vanished during the period 1975-2020 due to reduction in suitable habitats . Although its distribution throughout Europe is also highly fragmented, the northern dragonhead was classified as Least Concern in the most recent version of the IUCN Red List of Threatened Species (Ericsson et al., 2011). In Norway, northern dragonhead occurs exclusively in the southeastern part of the country. The most common habitats in Norway are dry, calcareous meadows, calcareous rocky outcrops along roads and railways, and extensively managed agricultural lands (Faegri & Danielsen, 1996; Norwegian Directorate for Nature Management, 2010).
It grows from a rhizome with limited vegetative branching, resulting in small clones of about 30-cm-tall stems with multi-flowered racemes. The peak season of the 2-to 2.5-cm-long, purplish-blue flowers is mid-June. Main pollinators are insects with long tongues, such as bumblebees (Apidae). As for most members of the mint family (Lamiaceae), however, northern dragonhead may also be selfcompatible (Milberg & Bertilsson, 1997). The species' generation time is approximately 15 years . An average fruit set rate of 0.27 has been reported for 20 populations in Sweden (Milberg & Bertilsson, 1997). Seeds are dry nutlets lacking modifications for long-distance dispersal and, hence, most likely are passively dispersed (Kyrkjeeide et al., 2020;Solstad et al., 2021).
Our selection of northern dragonhead as the target species for this study was based on four aspects: (1) The plant's biology as insectpollinated and mainly outcrossing (Milberg & Bertilsson, 1997). This is to avoid the complicating aspects of wind-pollination, polyploidy, and extensive selfing. (2) The abundance availability of specimens through time and space in Norwegian herbaria. (3) The availability of a 96 SNP microfluidic array for northern dragonhead (Kleven et al., 2019), which had already been used for analyzing contemporary samples from Norway (Kyrkjeeide et al., 2020). And, (4) the species' high level of priority in Norway, being threatened by habitat loss. Kyrkjeeide et al. (2020) added a genomic level to the monitoring regime of northern dragonhead in Norway. They found that only two of the revealed genetic groups were covered by the demographic monitoring and conservation efforts at the time.
An exploratory aspect of this study includes testing the performance of a microfluidic array for SNP genotyping, which has been developed based on genomic data from modern specimens of northern dragonhead, on up to 200 years old herbarium specimens of the species. Using this approach, we wanted to study the genetic diversity in northern dragonhead both back in time (i.e., prior to 1950) and across space (mainly in Norway but also beyond). More specifically, we wanted to assess whether the overall and regional genetic diversity of northern dragonhead in Norway has changed. Performing environmental niche modelling (ENM) on occurrence records of northern dragonhead in Norway, we wanted to reveal which abiotic factors may be limiting its distribution and whether there are areas in which northern dragonhead do not occur today that may represent suitable habitats.

| Materials
For the present study, we sampled 130 herbarium specimens of northern dragonhead collected between 1820 and 2008 (Table S1).
The majority of specimens originate from Norway, but we also included 25 specimens from Sweden, 12 from Russia, four from Ukraine, and two from each of the countries Belarus, Switzerland, F I G U R E 1 Map of Europe displaying the origins of all 130 included herbarium specimens of Dracocephalum ruyschiana. Each point represents a specimen colored according to the country of origin (see inset legend). NOR, Norway; SWE, Sweden; CHE, Switzerland; FRA, France; BLT, Belarus; UKR, Ukraine; RUS, Russia. and France (Figure 1). We additionally included published SNP data from 355 contemporary Norwegian samples of northern dragonhead from 43 different sites (Kyrkjeeide et al., 2022).

| DNA extraction and genotyping
All manipulations of the herbarium specimens postsampling were performed within the NTNU University Museum's dedicated, positively pressurized, paleo-genomics laboratory. About 0.5 cm 2 of leaf material was removed from each herbarium specimen using clean forceps and placed directly into a microfuge tube. The leaf material was pulverized with two Qiagen 3 mm tungsten carbide beads in a Qiagen TissueLyser LT at 50 Hz for 2 min. We extracted DNA using the DNeasy® Plant Mini Kit (Qiagen) with modifications from the manufacturer's instructions as described by Martin et al. (2014). We incubated the samples for 15 min at 37°C prior to spinning during the elution step. All extractions were performed using UV-sterilized equipment, and blank samples were always included to monitor for contamination. We measured DNA yield for 116 of the 130 herbarium samples using the Qubit dsDNA BR Assay Kit (Thermo Fisher Scientific, MA, USA), following the manufacturer's protocol. The DNA integrity was evaluated for the same 116 samples through agarose gel electrophoresis. The brightest band on the gel was regarded as an approximation for the sample's mean DNA fragment length. All samples were genotyped using the 96 × 96 SNP array developed by Kleven et al. (2019). The samples were genotyped on a Fluidigm EP1 instrument (Fluidigm Corporation, San Francisco, USA) according to the manufacturer's protocol and scored using the Fluidigm SNP genotyping analysis software v. 4.5.1 (https://www.fluid igm.com/ software). Positive and negative controls were included.
We excluded SNPs with more than 10% missing data across all herbarium samples. For each genotyped sample, we subsequently calculated the call rate (CR) and the proportion of successfully genotyped loci. The relationships between specimen collection year, CR, and DNA concentration were estimated using the Pearson's correlation test implemented in R (R Development Core Team, 2021).

| Genetic structure and diversity
The contemporary samples from Kyrkjeeide et al. (2022) are herein referred to as the modern group, whereas the herbarium samples are referred to as the historical group. We established two main datasets, NOR and GLOB, which we analyzed in their entirety or as subsets for assessing genetic structure and diversity. The NOR dataset included SNP data from Norwegian samples only, both historical and modern, and the former only from herbarium specimens collected prior to 1950 (N = 76). Due to incomplete overlap between historical and modern sampling sites (some even representing extinct populations), which complicated direct temporal comparisons of populations, we grouped samples into the geographical regions that correspond to the distinct genetic groups discovered by Kyrkjeeide et al. (2020).
Regions containing a sole historical sample were included only in our analyses of genetic structure. The GLOB dataset included SNP data from all 130 genotyped herbarium specimens, independent of collection year (83 Norwegian and 47 extra-Norwegian; Figure 1). The NOR and GLOB SNP datasets used in the present study are publicly available at DRYAD (https://doi.org/10.5061/dryad.c59zw3r8g).
The genetic structure within northern dragonhead was assessed using ParallelStructure v. 2.3.4 (Besnier & Glover, 2013;Pritchard et al., 2000), an R-based implementation of the common STRUCTURE algorithm, on XSEDE at CIPRES Science Gateway v. 3.1 (Miller et al., 2011). For the NOR dataset, we tested K from 1 to 40, with 20 replicates for each value of K. For the GLOB dataset, we per- Calculation of the optimal number of clusters and visualization of the results were obtained using find.cluster (R package adegenet v. 1.3-1; Jombart, 2008) and StructureSelector (Li & Liu, 2018), which implements the Puechmaille method of optimization (Puechmaille, 2016) in addition to calculating Ln Pr(X|K) and ΔK (Evanno et al., 2005).
For both the NOR and GLOB datasets, genetic structure and differentiation was evaluated through the use of discriminant analyses of principal components (DAPC; Jombart et al., 2010). Groups were defined a priori according to geography; by region for NOR and by the country for GLOB. For NOR, we also tested a priori grouping by historical versus modern samples. We additionally calculated the fixation index (F ST ) for NOR using the R function stamppFst (package StAMPP; Pembleton et al., 2013), after converting our data into a genlight object. Pairwise F ST values were calculated between the overall historical and modern groups, and between regional areas within the historical and modern groups, respectively, in addition to between the two age groups within the same regional areas. The 95% confidence interval was estimated using 1000 bootstraps.
For NOR, after transforming our data to a genind object, we calculated observed (H O ) and expected heterozygosity (H E ), and the inbreeding coefficient (F IS ) using the basic.stats function in R. We used the R summary function (adegenet package) to obtain the number of alleles, isPoly for the number of polymorphic loci, and private_alleles (R package poppr; Kamvar et al., 2014) for the number of private alleles. All measures of genetic diversity were calculated for the historical and modern regions, and averaged across all regions within each age group. In order to evaluate the effect of the sample size on the reported values, we subsampled modern regions down to equal sample size as their equivalent historical region. The subsampling was done randomly in 10 replicates using the R function sample, and diversity measures were recalculated in each run. Subsequently, we averaged across all replicates and calculated the standard deviation (SD). To visually explore the potential change in genetic diversity over time, we plotted the estimated H E of historical and modern regions against the oldest and youngest collection year, respectively.
The subsample averaged H E (and SD) was used for regions with uneven historical and modern sample size. We additionally plotted the individual proportion of heterozygous loci (PHt) against the sample collection year within each region.
For GLOB, genetic diversity was estimated as the proportion of heterozygous loci per individual (PHt). We calculated the number of polymorphic loci for the individual countries utilizing the adegenet function isPoly. To evaluate whether the GLOB genetic diversity was affected by distance from the SNP array's source population (Norway), we estimated the Pearson's correlation between PHt and the sample localities' distance from Oslo (59.9138 N, 10.7387E), an approximate centre of the SNP array's source populations.

| Environmental niche modeling
Species occurrence records for preserved specimens of Norwegian northern dragonhead were downloaded from the GBIF (10.15468/ dl.748g3v, accessed via GBIF.org on 2021-03-13). We added coordinates for GBIF-IDs lacking this information according to locality information and its precision (Table S2). Prior to analyses, we removed occurrence records that were clearly originating from a garden or otherwise represented a cultivar (Table S3).
Modeling the species' distribution was based on a final dataset of 4092 species occurrence records. The environmental niche modelling of northern dragonhead was based on three variables: mean temperature of the warmest quarter (i.e., mean summer temperature; MST), mean annual precipitation (MAP), and precipitation seasonality (coefficient of variance of monthly precipitation; PS). These are the same variables that, according to Speed and Austrheim (2017), represent the majority of uncorrelated variation in the total bioclimatic space of Norway. The climatic data were downloaded from WorldClim (Fick & Hijmans, 2017) at 1-km resolution (Figures S1-S12). We created background data by sampling 1000 random occurrence points across Norway, weighted by the distribution of occurrence data of all vascular plants in Norway. The sdm R package (Naimi & Araújo, 2016) was used to run several different distribution models: generalized lin-

| SNP genotyping performance
For the herbarium specimens, the DNA stock concentration varied from 1.17 to 41.30 ng/μl, with a mean of 16.11 ng/μl. The DNA stock concentration and specimen's collection year were positively correlated with a Pearson's correlation coefficient of R = 0.35 (p < .001; Figure 2a). Due to a technical error, four SNPs had missing data for more than 10% of all genotyped samples. We removed these prior to analyses. Across all historical samples and the remaining 92 SNPs screened, the mean call rate (CR) was 99.71%, ranging from 95.65% to 100%. In comparison, the modern samples had a mean CR of 99.89%, varying from 94.57% to 100%. When separating historical Norwegian versus extra-Norwegian samples, the mean CR was 99.91% and 99.35%, respectively. Considering only historical specimens, there was no significant correlation between CR and collection year (Pearson's correlation coefficient; R = 0.04 p = .68;

| Genetic structure and diversity -Norwegian scale
For the STRUCTURE analysis conducted on the NOR dataset (N = 431), the optimal number of genetic clusters varied depending on the applied optimization method. The mean log posterior ln P(K) was found to continuously increase with increasing K and reached the highest value for K = 24 ( Figure S2). We found the highest value of ΔK for K = 2, although ln P(K) was low at K = 1. Highest MedMed K, MedMean K, and MaxMean K were observed for K = 7, whilst MaxMed K was highest for K = 8. Using find.cluster, the lowest BIC value was found between K = 5 and K = 8 ( Figure S3). Under the most frequently inferred number of clusters (K = 7), when sorting samples according to predefined regions, six of the clusters largely corresponded to the regional areas: Hedmark, Oslofjorden (east and west), Randsfjorden, Tyrifjorden, and Valdres-Gudbrandsdalen ( Figure 3a). The other regional areas that appeared admixed for all Ks (i.e., Agder, Drammensfjorden, Hemsedal, and Ytre Oslofjorden; Figure S4) were excluded from downstream analyses because each contained only a single sample. When dividing our STRUCTURE results into historical and modern groups, most of the regional areas displayed similar genetic structures through time (Figure 3b). The greatest temporal change in ancestry proportions was observed for Randsfjorden, whereas the least change through time was observed for Oslofjorden.
The DAPC results corroborated the separation of Oslofjorden and Tyrifjorden, respectively, from the remaining regions, when a priori grouping our samples according to the regional areas ( Figure 4).
The first and second DA explained 53.4% and 24.7% of the genetic variation, respectively. A priori grouping of the specimens by age (historical vs. modern) for the DAPC analysis resulted in overlapping density curves along the first axis ( Figure S5). Also, in terms of F-statistics, we observed larger genetic divergence across space than through time. The overall F ST value between the historical and modern groups indicated a very low overall level of temporal genetic divergence (i.e., 0.003 with 95% CI from 0.002 to 0.004). Pairwise comparisons of regions (in both time and space) yielded generally low F ST values, but all 95% confidence intervals were above zero (Table S4). The largest temporal genetic divergence was observed for Randsfjorden (F ST = 0.033) and Tyrifjorden (F ST = 0.028), which both also displayed a decline in spatial divergence over time. This decline in F ST was only significant for Randsfjorden (Table 1) Table 2). The number of alleles and polymorphic loci was largest within the modern group ( Table 2). No private alleles were found for the historical versus modern group, and no specific region contained private alleles compared to the other regions within the same age group. We did, however, observe private alleles when comparing the historical and modern samples within single regions ( Table 2). The highest amounts of private alleles were found in modern Randsfjorden (25) and Tyrifjorden (21) compared with their respective historical regions, likely a result of uneven sample sizes F I G U R E 2 Correlation plots for DNA quality (as stock concentration and mean fragment size), call rate (CR; the proportion of successfully genotyped SNPs per sample), and age of the Dracocephalum ruyschiana samples (given as collection year). The orange lines represent the average overall historic samples and the orange zones the 95% confidence interval. Each symbol represents an individual sample, the shape of its geographical origin, and the color of its mean DNA fragment size (bp) based on gel electrophoresis (white = no data).

| Genetic structure and diversity -Global scale
Results from the STRUCTURE analysis on the full GLOB dataset varied with regard to the number of optimal genetic clusters depending on the applied optimization method ( Figure S6). The mean log posterior, ln P(K), increased until K = 3 and made a drop before increasing to its maximum at K = 7. We found the highest value of ΔK under K = 2, although ln P(K) was low at K = 1. We found the highest MedMed K, MedMean K, MaxMean K, and MaxMed K for K = 4. Using find.cluster, we observed the lowest BIC value between K = 2 and K = 6 ( Figure S7). At K = 2, Norwegian samples separated from the remaining European samples ( Figure S8). At K = 4, Swedish samples formed their own group while French and Swiss samples displayed mixed ancestry from the Norwegian and Swedish clusters ( Figures S8-S8). Further increasing K led to a higher degree of admixture, mainly within Norway, but also to some degree within Sweden, Switzerland, and France ( Figure S8). Belarus, Russia, and Ukraine, on the other hand, consistently formed a single cluster. For the other two STRUCTURE analyses, excluding Norwegian samples (Figures S10-S11) and balancing sampling across countries (Figures S12-S13), we observed that ln P(K) increased until K = 4 and K = 3, respectively. For increasing values of K, the ln P(K) continued to decrease. For both these analyses, the highest value of ΔK was found under K = 2, whereas MedMed K, MedMean K, MaxMean K, and MaxMed were highest at K = 3.
When a priori grouping the herbarium specimens by geography (i.e., by country), the first and second DA explained 59.3% and 40.7% of their total genetic variation, respectively ( Figure 6). The DAPC analysis separated the Norwegian population from the remaining Eurasian countries. We also found the highest genetic diversity within the Norwegian samples, measured as individual proportions of heterozygosity (PHt) and the number of polymorphic loci ( Figure S14a, b). The individual PHt decreased significantly with increasing distance from Norway (R = −0.49, p = 5.69 × 10 −9 ; Figure S14c).

| Environmental niche modeling
Across all the replicated environmental niche models, the mean AUC was 0.95 (SD ± 0.03). The relative variable importance was highest for mean summer temperature (MST; 0.40, SD ± 0.02), followed by mean annual precipitation (MAP; 0.31, SD ± 0.02), and lowest for precipitation seasonality (PS; 0.15, SD ± 0.01; Figure 7a). Based on the response curves, climate suitability for northern dragonhead increased with higher temperatures, MST > 10°C, and decreased with increased precipitation, MAP > 500 mm (Figure 7b). After averaging over all models, the model predicted the greatest niche suitability in southeastern Norway (Figure 7c). Potentially suitable, but unoccupied, niche space was predicted around Trysil and in the lowland valleys east of the Trondheimsfjord, among other areas in western and northern Norway (Figure 7c).

| DISCUSS ION
Maintenance of genetic diversity is a central aim of species conservation, given its positive role in a species' performance and survival in a changing environment (Lande & Shannon, 1996). In this study, we have assessed changes in genetic structure and diversity across space and through time in northern dragonhead, a charismatic flowering plant that has experienced a drastic population decline and habitat loss in Europe. We have added a temporal level to the monitoring of northern dragonhead in Norway using an SNP array technology on herbarium specimens. To identify which abiotic factors may limit its distribution and whether there are additional areas with suitable habitats, we have used sample metadata and observational occurrence records to model the species' environmental niche.

| SNP genotyping performance through time
All the included herbarium specimens of northern dragonhead yielded DNA of a quality suitable for SNP genotyping. Even though the DNA stock concentration decreased with specimen age (Figure 2a), the negative correlation was weaker than expected.
Previous time-series studies of herbarium samples have shown that both molecular weight (DNA fragment length) and stock concentration decreased with time since collection (see Raxworthy & Smith, 2021, and references therein). The rate of decrease in molecular weight and DNA concentration apparently depends on the samples' history, such as the way it was collected and preserved, and the subsequent storage conditions. In addition, DNA concentration appears to vary among different parts of the specimen, tissue types, and preservation techniques. Indeed, for herbarium specimens, most of the DNA damage appears to occur soon after sampling (i.e., during specimen preparation; Staats et al., 2011). The best practice for preserving plant DNA is assumed to be rapid desiccation under moderate temperatures.  Note: The fixation index values (F ST ) represent pairwise comparisons of either different regional areas in historical times (yellow, lower triangle) or modern times (orange, upper triangle), or between the modern and historical groups within the same regional area (white, diagonal). The F ST values that have changed significantly are marked with an asterisk. Only one asterisk means that the value is lower than that of the other age group, and two asterisks means it is higher than that of the other age group. The 95% confidence intervals are displayed in Table S4.  Bendiksby et al., 2011). Moreover, the fact that our study object occurs in the temperate zone, rather than in the tropics, implies that the specimens studied have been living in a less harsh climate (i.e., moderate temperatures) with better facilities for rapid desiccation.

F I G U R E 4 DAPC results for our Norwegian
Hence, although the microfluidic SNP array approach was highly successful for northern dragonhead, this may not be the case for historical specimens of species that, for biological reasons, experience faster DNA degradation, or that cannot be desiccated rapidly under moderate temperatures.

| SNP genotyping performance across space
For the global dataset (GLOB), our results demonstrate the effect of SNP ascertainment bias (i.e., the selection of loci from an unrepresentative sample of individuals), which shows a systematic deviation from theoretical expectations (Geibel et al., 2021).
Since the SNP array was designed based on highly polymorphic SNPs from Norwegian northern dragonhead samples, the allele frequencies were expected to be lower in populations outside Norway. This is apparent from our global measures of genetic diversity, which decrease significantly with increasing geographical distance from Norway ( Figure S14c). On the other hand, ascertainment bias is apparently less likely to affect the assignment of individuals to separate populations (Lachance & Tishkoff, 2013).
As such, our results indicate that the Norwegian samples are genetically distinct from the examined materials originating from elsewhere in Eurasia (Figures 6, S9). To determine the degree to which they are distinct cannot, however, be estimated based on our current SNP data.

| Microfluidic SNP array optimization
The critical step for obtaining informative SNP data lies in the selection of SNP markers and the development of the SNP array itself.
Since genetic diversity is often unevenly distributed across space and through time, the SNP data will be biased towards variants present in the samples from which the selected SNPs originate (Geibel et al., 2021). For example, genetic diversity only present in the past will not be recovered by an array designed based on modern material alone. To reduce ascertainment bias and to obtain more precise genetic estimates for a spatiotemporal study, the array of SNPs should

| SNP array + herbarium = cost-and risk savings
Apart from successfully genotyping historical herbarium specimens, SNP genotyping with microfluidic arrays also offers a cost-and time-efficient method for generating genomic datasets for many samples (von Thaden et al., 2017). This is particularly the case when genomic data for SNP selection is already available (e.g., genome skims or RAD/GBS data). Prior to loading the DNA onto the array, no library preparation is required, and large numbers of samples can be processed simultaneously. Furthermore, there is no need for extensive bioinformatic skills, the raw data require less storage space, and the computational time for filtering and processing the data is comparably short. Hence, for processing many samples for genetic monitoring purposes, including historical ones with variously degraded DNA for temporal monitoring, microfluidic SNP genotyping appears to be a promising method of choice due to reduced overall cost and labour as compared to other currently available methods.
Working with historical specimens provides several key benefits compared to using contemporary material alone. As demonstrated by our study, incorporating historical specimens, which could even include extinct populations, enables the assessment of genetic We did record a small decrease in the overall inbreeding coefficient over time, from F IS = 0.107 in the historical to F IS = 0.062 in the modern group. A decrease in F IS (but still F IS >0) could be an indication of overall less effects of genetic drift or a higher degree of outcrossing compared with historical times. It should, however, be mentioned that the standard deviation for the obtained overall F IS was relatively large for the historical group (SD ± 0.056).
Moreover, our current data may not be suitable for robustly infer- Slovakia (Dostálek et al., 2010).
A seemingly unchanged distributional range in Norway and limited dispersal may also have contributed to the observed temporal genetic "stasis" within northern dragonhead. Despite a reduction in suitable habitats over the last 150 years, observational data indicate that the overall distributional range of northern dragonhead in Norway has remained largely intact, and that the decline has been mainly local rather than regional (Norwegian Directorate for Nature Management, 2010: Figure 6). Given its pollination syndrome (insect pollination; Milberg & Bertilsson, 1997) and relatively large seeds, the northern dragonhead is primarily an outcrossing species with presumably poor abilities for long-distance dispersal. Additionally, the landscape topology of Norway, corresponding well with the predefined regional areas used herein (adopted from Kyrkjeeide et al., 2020), likely limits dispersal between regions naturally.
It should be mentioned that isolation by distance (IBD) was shown to be present in modern samples of northern dragonhead in Norway (Kyrkjeeide et al., 2020: Figure 2). Their Mantel test revealed a positive correlation between genetic distance and geographical distance (R = 0.56, p = .001). The analysis software STRUCTURE, which we have used herein, assumes that markers are not linked and that populations are panmictic (Pritchard et al., 2000). Hence, our STRUCTURE results should be interpreted with caution, as IBD violates the assumption of freely distributed genotypes. In our study, however, also the DAPC results support that genetic variation within historical and modern northern dragonhead is better explained by divergence across space than divergence through time. The DAPC analysis software is a model-free method based on K-means clustering of genetic distance and IBD does not violate its assumptions (Jombart et al., 2010). 4.2.2 | Minor temporal genetic change at regional level At the regional scale, the temporal genetic changes were also small. The direction of change, however, varied between regions.
For four of the regions (Gudbrandsdalen, Oslofjorden, Tyrifjorden, and Valdres), the inbreeding coefficient decreased over time ( Table 2). There was, however, still an excess of homozygosity rela-  and <~800 mm of mean annual precipitation (Figure 7b). These findings are in line with the early assumption by Sterner (1922) that the distribution of northern dragonhead is limited by low summer temperatures. In its current southeastern distribution in Norway, northern dragonhead is further restricted to areas of dry, calcareous meadows or steep, rough land like ledges along roads, in addition to extensively managed agricultural lands (Faegri & Danielsen, 1996). Further, east of its present distribution, the valleys are dominated by noncalcareous soils and bedrock not suitable for northern dragonhead (Faegri & Danielsen, 1996). However, our ENM results suggest areas representing potentially suitable climatic niche space for northern dragonhead in Trøndelag (central Norway), the inner parts of the fjords in the western part of the country, and in northeastern Norway (Figure 7c, deep red).
The latter area may seem unlikely given the cold and long winters above the Arctic circle at approximately 70 degrees north. Notably, this area was suggested as suitable also for Carex jemtlandica (see Nygaard et al., 2021), which also has a mainly southeastern

| CON CLUS ION
With this study, we demonstrate that, with the appropriate design procedures, the microfluidic SNP array technology is promising for genotyping old herbarium specimens; an invaluable source of information from the past. As expected, the SNP array picked up less genetic variability in the extra-Norwegian specimens, likely due to both genetic divergence and the fact that the array was developed based on modern Norwegian samples alone. Our temporal genomic analyses of northern dragonhead in Norway show no signs of any severe reduction in population size in any of the studied regions. This may seem like good news, which indeed it might be if it is so that the populations have remained large enough to withstand the effect of genetic drift and inbreeding. The same results may, however, be due to a time lag in the response caused by the relatively long generation time of northern dragonhead. It is tempting to speculate whether our results could also be reflecting the ongoing climate change; increasing temperatures and less precipitation could potentially lead to an increase in connectivity and gene flow between neighboring populations and an expansion of the limits of currently suitable habitats. Regardless, the regional areas studied are genetically divergent across space, both from each other and clearly so from populations outside of Norway, rendering continued protection of the species and its regional genetic variation in Norway relevant. Our ENM results suggest that northern dragonhead has not yet reached its potential distribution in Norway. With the future inclusion of additional parameters (e.g., pH), ENM should prove useful for guiding management authorities in translocation for conservation initiatives.

ACK N OWLED G M ENTS
We acknowledge Hans Stenøien and Vibekke Vange for taking part in an early brainstorming of ideas, of which some have become part of this publication. MN thanks Vanessa Bieker for training in the ancient DNA lab facility at NTNU University museum and for assisting with the first DNA extractions. We specially thank Mikhail P. Zhurbenko for providing Russian specimens. We thank the curators at O, TRH, and UPS for letting us borrow and make use of their specimens. We are grateful to all the collectors of the specimens we have studied. We also thank colleagues at NINA for letting us reuse their published SNP data from modern specimens of Dracocephalum ruyschiana (Kyrkjeeide et al., 2020). Seed money from the NTNU University Museum, Norwegian University of Science and Technology, Trondheim, Norway, funded this project.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest and no relevant financial or nonfinancial interests to disclose.

DATA AVA I L A B I L I T Y S TAT E M E N T
The herbarium specimens used for the newly generated data are in the following public collections: LECB, O, TRH, and UPS, and voucher information provided in Table S1. The species occurrence records used herein are available from GBIF (see the link provided in chapter 2.4), and Tables S2 and S3 list the GBIF-IDs lacking coordinates (i.e., our approximations based on specimen metadata) and those that were considered spurious, respectively. Both our newly produced SNP data and those produced by Kyrkjeeide et al. (2020)